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ABSTRACT 



Context. The asteroseismology of rapidly rotating pulsating stars is hindered by our poor knowledge of the effect of the rotation on 
the oscillation properties. 

Aims. Here we present an asymptotic analysis of high-frequency acoustic modes in rapidly rotating stars. 

Methods. We study the Hamiltonian dynamics of acoustic rays in uniformly rotating polytropic stars and show that the phase space 
structure has a mixed character, regions of chaotic trajectories coexisting with stable structures like island chains or invariant tori. In 
order to interpret the ray dynamics in terms of acoustic mode properties, we then use tools and concepts developed in the context of 
quantum physics. 

Results. Accordingly, the high-frequency acoustic spectrum is a superposition of frequency subsets associated with dynamically 
independent phase space regions. The sub-spectra associated with stable structures are regular and can be modelled through EBK 
quantization methods while those associated with chaotic regions are irregular but with generic statistical properties. The results of 
this asymptotic analysis are successfully confronted with the properties of numerically computed high-frequency acoustic modes. The 
implications for the asteroseismology of rapidly rotating stars are discussed. 
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' 1. Introduction 

■ Interpreting the oscillation spectra of rapidly rotating stars is a 
' long standing unsolved problem of asteroseismology. It so far 
prevented to directly probe the interior of stars in large frac- 
tions of the HR diagram, mostly in the range of intermediate 
and massive stars. Important progresses are nevertheless ex- 
pected from the current spatial seismology missions (in partic- 
ular COROT and KEPLER) as well as from recent modelling 
efforts on the effect of rotation on stellar oscillations. On the 
modelling side, new numerical codes have been able to accu- 
rately compute oscillation frequencies in significantly centrifu- 
' gaily distorted po lytropic models of stars (iLignieres et al.Ll2006t 
iRee se et aL. 2006) and are being extended to more realistic mod- 
els GReese et a l., 2009 ). In particular the previous calculations 
based on pertu rbative methods were shown to be inadequat e 
for these stars dReese et all 120061: iLovekin & Deupreel l2008l) . 
Nevertheless, interpreting the data requires an understanding of 
the mode properties that goes far beyond the accurate compu- 
tation of frequencies. Indeed, the necessary identification of the 
observed frequencies with theoretical frequencies is a largely un- 
derconstrained problem which requires a priori information on 
the spectrum to be successful. The knowledge of the structure of 
the frequency spectrum is of primary importance in this context. 
For slowly rotating pulsating stars, this structure is character- 
ized by regular frequency patterns which can be analytically de- 
rived from an asymptotic theory of the high-frequency acoustic 
modes. 
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Until recently, the asymptotic structure of the frequency 
spectrum of rapidly rotating stars was completely unknown. Our 
first calculations of low degree (f - 0-7) and low order 
(n - 1-10) acoustic axisymmetric modes in centrifugally dis- 
torted polytropic stars (ILignieres et al.L l2006h have revealed the 
presence of regular frequency patterns similar albeit different 
from those of spherically symmetric stars. This was confirmed 
using more realistic calculations including the Coriolis force and 
was als o extended t o non -axisymmetric and higher-frequency 
modes (i Reese et"an, l2008h . The analogy with the non-rotating 
case suggests an asymptotic analysis could model these empiri- 
cal regular patterns. 

The asymptotic analysis presented in this paper is based 
on the acoustic ray dynamics. This approach can be viewed 
as a natural extensio n of the asymptotic analys i s developed 
for non-r otating stars (Vandakurov', '1967"; 'Tassoul', ' 198(1, 1 19901; 
iDeubner & Gough, 1984; Roxburgh & Vorontsov, 2000). In this 
case, spherical symmetry enables to reduces the initial 3D 
boundary value problem to a ID boundary value problem in the 
radial direction. Asymptotic solutions of this ID boundary value 
problem can then be obtained using a short-wavelength approx- 
imation which consists in looking for wave-like solutions un- 
der the assumption that their wavelength is much smaller than 
the typical lengthscale of the background medium. As rotation 
breaks the spherical symmetry, the eigenmodes are no longer 
separable in the latitudinal and radial directions and the 3D 
boundary value problem of acoustic modes in rapidly rotating 
stars can not be reduced to a ID boundary value problem. Thus, 
the short-wavelength approximation is directly applied to the 3D 
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equations governing linear adiabatic stellar perturbations. This 
leads to an acoustic ray model which describes the propagation 
of locally plane waves. Since we are concerned by oscillation 
modes, the main issue of an asymptotic analysis based on ray dy- 
namics is to construct standing waves solutions from the short- 
wavelength propagating waves described by the acoustic rays. 

The short-wavelength approximation of wave equations is 
standard in physics, best known examples being the geometrical 
optics limit of electromagnetic waves or the classical limit of 
quantum mechanics. We shall see that similarly to these cases 
the acoustic rays in stars can be described as trajectories of a 
particle in the framework of classical Hamiltonian mechanics. 
As is well known in quantum physics (Gutz wilier,, .1990; Ott, 
fT993h . the issue of construction modes from the ray dynamics 
depends crucially on the nature of this Hamiltonian motion. 

Indeed, Hamiltonian systems can have one of two very dif- 
ferent behaviors. If there are enough constants of motion, the 
dynamics is integrable, and trajectories organize themselves in 
families on well-defined phase space structures. In contrast, 
chaotic systems display exponential divergence of nearby trajec- 
tories, and a typical orbit is ergodic in phase space. The modes 
constructed from these different dynamics are markedly differ- 
ent. For integrable systems, the eigenfrequency spectrum can be 
described by a function of integers, being the number of 
degrees of freedom of the system. In contrast, the spectrum of 
chaotic systems shows no such regularities. However, the spec- 
trum possesses generic statistical properties which can be pre- 
dicted. In the last thirty years, the field called quantum chaos 
has developed concepts and methods to relate non-integrable ray 
dynamics to properties of the associated quantum systems (and 
more generally wave systems). 

We shall see that the acoustic ray dynamics in rotating stars 
undergoes a transition from an integrable system at zero rotation 
to a mixed system, that is a system with a phase space containing 
integrable and chaotic regions, at high rotation. The acoustic ray 
dynamics of rotating stars being non-integrable, we are led to 
use quantum chaos theory to predict the asymptotic properties 
of acoustic modes. 

In the following, we shall describe in detail the ray dynam- 
ics, the predictions on the modes properties and their valida- 
tion through a comparison with numerically computed acous- 
tic modes. But, before going into the technical details of this 
analysis, we would like to summarize our results, giving empha- 
sis to those of practical interest for mode identification. These 
results have been obtained for a sequence of uniformly rotat- 
ing polytropic models, but we expect them to be qualitatively 
correct for a wider range of models. At high rotation rates, the 
frequency spectrum can be generically described as the super- 
position of an irregular frequency subset and of various regu- 
lar frequency subsets each showing specific patterns. This spec- 
trum structure is significantly more complex that in the spherical 
case where all acoustic frequencies fo llow the s ame o rganization 
such as given by Tassoul's formula (iTassoull Il980b . However, 
in the observed spectrum, only two frequency subsets are ex- 
pected to be dominant. One subset (the subset of island modes) 
shows regular patterns similar albeit different to those found 
in no n-rotat ing stars (it corres ponds to t he modes subset stud- 
ied bv lLignieres et al.l (120061) : lReese et alj (12008 )). These modes 
are associated with rays whose dynamics is near-integrable and 
consequently asymptotic formulas describing their regular pat- 
terns can be derived. The second frequency subset (the subset 
of chaotic modes) is associated with chaotic rays. Although it 
does not follow a regular pattern, specific statistical properties 
of this frequency subset can be predicted. Besides, the asymp- 



totic analysis provides an estimate of the proportion of mode in 
each subset. The transition from the non-rotating case occurs as 
follows: At moderate rotation, the regular subset of island modes 
is superposed to another regular subset (the subset of whispering 
gallery modes) which is a direct continuation of the non-rotating 
spectrum. At this stage, chaotic modes are rare and difficult to 
observe. As rotation increases, the number of chaotic modes in- 
creases while whispering gallery modes become less and less 
visible. Obviously, such a priori informations on the frequency 
spectrum should be of practical interest for the mode identifi- 
cation. One should however keep in mind that the asymptotic 
analysis is not supposed to be accurate for the lowest frequency 
p-modes. Although a detailed study of the asymptotic analysis 
limit of validity has not been performed yet, numerical results 
indicate that the regul ar patterns are quite accurate down to 5* 
radial order (see Ligni eres et al.L 2006; Reese, 2007). At lower 
radial orders, the asymptotic mode classification in different sub- 
sets could still be applicable. 

The paper is organized as follows: The equations govern- 
ing the star model, the small perturbations about this model, 
the short-wavelength approximation of these perturbations and 
the ray model for progressive acoustic waves are presented in 
Section 2. A detailed numerical study of the acoustic ray dynam- 
ics has been conducted for uniformly rotating polytropic models 
of stars. 

The results are analyzed in Section 3 using Poincare Surface 
of Section, a standard tool to visualize the phase space structure. 
It shows that, as rotation increases, the dynamics undergoes a 
transition from integrability (at zero rotation) to a mixed state 
where parts of the phase space display integrable behaviour and 
while other parts are chaotic. 

We then relate the acoustic ray dynamics to the asymptotic 
properties of the acoustic modes (Section 4). We first summarize 
the results obtained in the context of quantum physics distin- 
guishing the cases where the Hamiltonian system is integrable, 
fully chaotic or of mixed nature. In accordance with this theory, 
we show that the asymptotic acoustic spectrum of the uniformly 
rotating polytropic models of stars is a superposition of regu- 
lar frequency patterns and irregular frequency subsets respec- 
tively associated with near-integrable and chaotic phase space 
regions. The averaged number of modes associated with each 
phase space regions depends directly on its volume (in phase 
space). These predictions are then successfully confronted with 
the actual properties of high-frequency acoustic modes com- 
puted for a particular fast rotating stellar model 

In Section 5, after a critical discussion of the assumptions of 
the asymptotic analysis, we show how our results can be used 
for the mode identification and for the seismic studies of rapidly 
rotating stars. The conclusion is given in section 6. 

Th e present work complemen ts and extends a short recent 
paper dLignieres & Georgeotl l2008i) by giving all the necessary 
details of the analysis and by presenting new results concern- 
ing (i) the ray dynamics at different rotation rates and for non- 
vanishing values of the angular momentum projection onto the 
rotation axis L., (ii) the analysis of extra regular patterns visible 
for some specific values of rotation, (iii) the number of modes in 
each frequency subset predicted by the asymptotic analysis, (iv) 
the visibility of the chaotic modes. 

2. Formalism and numerical methods 

In this section we present the equations governing the star model 
(subsection 2.1), the small perturbations about this model (sub- 
section 2.2), the short-wavelength approximation of these per- 
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turbations (subsection 2.3) and the ray model for progressive 
acoustic waves (subsection 2.4). The numerical method used to 
compute the ray trajectories is described in subsection 2.5. The 
os cillation modes have b een computed with the code described 
in lLignieres et al.l (l2006l) and lReese et al. I (l2006l) . 

2.1. Polytropic model of rotating star 

The model is a self-gravitating uniformly rotating monatomic 
gas (F = 5/3) verifying a polytropic relation assumed to give a 
reasonably good approximation of the relation between the pres- 
sure and the density in the star (lHansenlil994i) : 
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where Pq is the pressure, po the density, K the polytropic con- 
stant, fj. the polytropic index, i/tq the gravitational potential, Q 
the rotation rate, w the distance to the rotation axis and G the 
gravitational constant. 

The uniform rotation ensures that the fluid is barotropic. A 
pseudo-enthalpy can then be introduced Iiq - J dPo/po = (1 -H 
yu)f o/po and the integration of the hydrostatic equation reads: 
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2... 2 



ho = he - (lAo - lAc) + ■:^n w 



(4) 



where the subscript "c" denotes the value in the center of the 
polytropic model. Equation (|4|i is then inserted into Poisson's 
equation to yield: 



Ai/Co = AnGpc 1 - 



!,. 2hr 



(5) 



Equation ^ is s olved numerically, u sing an iterative 
scheme, as described in lRieutord et alj (l2005h . 

Specifying the mass and the rotation rate of the star is not 
sufficient to determine the polytropic model in physical units. 
This requires to fix an ad ditional parameter, for example the 
stellar radius (see HansenI ( 1994.) for the no n-rotating case and 



IChristensen-Dalsgaard & ThompsonI (Il999h for a brief discus- 
sion on a suitable parameter choice for rotating stars). In the 
following, however, we only present dimensionless quantities 
which do not depend on the choice of this additional parameter. 

/ •,\l/2 

The rotation rate Q is compared to Q^: = [GMIRA the lim- 
iting rotation rate for which the centrifugal acceleration equals 
the gravity at the equator, M being the stellar mass and Rg the 
equatorial radius. The star flatness is e = 1 - Rp/Rg where Rp 
is the polar radius. The acoustic frequencies shall be expressed 

in terms of ojo - (GM/R^j^ ^ , the inverse of a dynamical time 
scale, or wi the lowest acoustic mode frequency of the stellar 
model considered. 

2.2. Perturbation equations and boundary conditions 

Time-harmonic small amplitude perturbations of the star model 
are studied under two basic assumptions. The first is to neglect 
the Coriolis force. This a natural assumption to study high- 
frequency acoustic modes since the oscillation time scale is 
asymptotically much shorter than the Corio lis force time scale 
l/(2f2). Moreover, complete calculations bv lReese et al. I (f2006h 
(see figure 6 of this paper) have shown that Coriolis force effects 



on the frequency are already very small for relatively low radial 
order (say n « 5). The second basic assumption is to neglect 
the viscosity and the non-adiabatic effects. This is a standard ap- 
proximation in the asymptotic analysis since these effects have 
little influence on the value of the frequency. Both assumptions 
shall have important consequences on the acoustic ray dynamics 
described below. Neglecting the Coriolis force ensures that the 
dynamics is symmetric with respect to the time reversal while 
the absence of diffusion processes makes the dynamics conser- 
vative. Finally, the Cowling approximation valid for high fre- 
quencies enables to neglect the perturbation of the gravitational 
potential. Under these assumptions, the linear equations govern- 
ing the evolution of small amplitude perturbations read: 



d,p + V ■ (pom) = 0, 

po<9,H ^-VP+pgo, 

d,P + u- VPo = cl (d,p + u ■ Vpo) , 



(6) 
(7) 
(8) 



where u, p and P, are respectively the Eulerian perturbations 
of velocity, density and pressure. The sound speed is Cj = 
^jTPolpa, r being the first adiabatic exponent of the gas, and 
i^o - Q ^ w^/2) is the effective gravity. 
PekerisI (119381) . because the pressure and the tempera- 
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= -V 
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ture of the stellar model is zero at the surface, the only condition 
to impose on the perturbations is to be regular everywhere. 

2.3. Tiie sliort-waveiengtli approximation of tiie perturbation 
equations 

The acoustic ray model results from a short-wavelength ap- 
proximation of the perturbation equations (|6]), d?), called 
the Wentzel-fCramers-Brillouin (WKB) approximation. Time- 
harmonic wave-Uke solutions of the form. 



= !R{A(x) exp[/0(x) - iujt]} = !R{4'(x) exp(-!wf)) 



(9) 



are sought under the assumption that their wavelength is much 
smaller than the typical length scale of the background medium. 
As discussed by lGoughl (Il993h . one expects a better approxima- 
tion if the starting equations (|6]l, d?), ^ are first reduced to a 
so-called normal form which avoids first order derivatives. This 
is done in Appendix A. 1 leading to: 
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(10) 



where 4* - P/a, P is the complex amplitude associated with the 
pressure perturbation P = ^R{P exp(-/dLtf)) and a is a function of 
the background star model given by Eq. ( IA.8b . The expressions 
of the cut-off frequency (li^ and the Brunt-Vaisala frequency A^o 
are given respectively by ( IA.5b and ( IA.3I ). Note that the two left- 
hand-side terms of equation ( fTOb account respectively for acous- 
tic and gravity waves. 

As detailed in Appendix A. 2, the WKB approximation is 
then applied to ( fTOl i. The dominant term of the expansion in 
powers of the ratio between the wavelength solution and the 
background typical lengthscale yields an equation governing the 
phase 'l'(x), the so-called eikonal equation. The amplitude A(x) 
is determined at the next order as a function of 0(x). Neglecting 
the gravity waves by taking the high frequency limit, the eikonal 
equation reads: 



2 2 , 2;2 



(11) 
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where k - VO is the wavevector Moreover, for a polytropic 
model of star and using the approximation oj » A^o valid for 
high frequency acoustic modes, the function a is proportional to 



1 /2 

Pg and the cut-off frequency cl>c is simplified to: 



-2 



2/10 



ji+l 



■ go. 



(12) 



In the range of high frequency acoustic modes, the w^, term is 
expected to be much smaller than oj in most parts of the star 
except near the surface where ojc diverges. Note also that, despite 
the chosen notation, a>l can be negative near the center. 

2.4. The acoustic ray model as a Hamiltonian system 

The acoustic ray model consists in finding solutions of the 
eikonal equation (fTTT i along some path called the ray path. This 
problem can be written in a Hamiltonian form where the so- 
lutions (x{t), kit)) are conjugate variables of the Hamiltonian 
and t, the parameter along the path, is a time-like variable. To 
derive the Hamiltonian equations from the eikonal equation, 
one can apply a pro cedure vah d for a general eikonal equation 
D{k,u,x) = (e.g. lOtl [T993I) . This leads to the Hamiltonian 
//' = w = Vc?k2 + 0)1 (e.g. iLighthilli Il978l) . Another useful 
Hamiltonian formulation can be readily obtained by considering 
a path normal to the wave front (^>{x) - const, (this method being 
also used to determine optical rays in isotropic media of variable 
index jBornl Il999l) . the quantity ( ^jl - (J^ I (jiP-) I c s playing the 
role of the medium index n{x)). The ray path is thus defined by 



dx _ k _ VO) 



(13) 



where s is the curvilinear coordinate along the ray. 
Differentiating ( fT3b and using ( fTTT l then leads to the following 
system of ODEs: 



dx ~ 

T ^ ^ 
dt 

dk 

dt 



2cl 



1 - 



(14) 



(15) 



(16) 



where we use the frequency-scaled wavevector k - kjio and the 
time-like variable f such that dt - Csdsjil - w^/w^)'^^. 

As W only depends on the spatial variable x, the second 
equation has the classical form of Newton's second law for the 
conservative force associated with the potential W (for a unit 
mass and a time variable t). This system can thus be written in a 
Hamiltonian form where 

H^ — + W{x) (17) 

is the Hamiltonian. According to the eikonal equation ( fTTT l. H 
is equal to zero and the dynamics is therefore fully determined 
by the potential well W, where the frequency w appears as a 
parameter As » Uc away from the near surface layers, the 
potential increase towards the star center is given by the sound 
speed increase. Close to the surface, the potential barrier is due to 
the sharp increase of Uc and provokes the wave back-reflection. 
While the location of the reflection surface depends on the fre- 
quency (x>, we shall see that the dynamics is not strongly depen- 
dent on (jj in the range of high-frequency acoustic modes here 



considered. This can be expected since, as Uc diverges towards 
the surface of the polytropic model of star, the location of w = w^. 
does not vary rapidly with to. 

The potential being symmetric with respect to the rotation 
axis of the star, the angular momentum projection on this axis 
L, - r sin 9k^ is a constant of motion, where - k ■ and is 
a unit vector in the azimuthal direction and [r, 6, <p] are the spher- 
ical coordinates. Consequently, the projection of the ray trajec- 
tory on the meridional plane rotating with the ray at an angular 
velocity d(p/dt = L~/(rsmd)^ is governed by the two-degree-of- 
freedom Hamiltonian; 



k 2 



(18) 



where kp = k-k^e^ is the wavevector projected onto the rotating 
meridional plane and Wr is the effective potential of the reduced 
Hamiltonian H,- which now also depends on as a parameter: 
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2.5. Numerical method for the ray dynamics 

The acoustic ray dynamics has been investigated by integrating 
numerically Eqs. ( fT4b and ( fTSl ) using a 5th order Runge-Kutta 
method. The step size of the integration is chosen automatically 
to keep the local error estimate smaller than a given tolerance. 
To what extent this control of the local error ensures that the 
numerical solution is close to the true solution depends on the 
stability of the problem. For chaotic trajectories, the numerical 
error tends to grow rapidly while for stable trajectories this error 
remains of the same order as the relative error. The rapid growth 
of numerical errors is one of the characteristics of chaotic dy- 
namics; however, this does not prevent to simulate such systems 
since for hyperbo lic systems the shadowing theorem (Anoso^ 
119671 ISauer et all 119971) ensures that an exact trajectory will 
remain close to the dynamics of each computed point for ar- 
bitrary times. Thus while a numerical trajectory diverges from 
the exact one, it nevertheless remains close to another exact tra- 
jectory, and therefore numerical errors do not prevent to obtain 
accurate phase space plots. We checked that the Poincare sur- 
face of section shown in the next section are not significantly 
modified by decreasing the maximum allowed local error We 
also checked the influence of the resolution of the background 
polytropic stellar model. Increasing this resolution from 62 to 
92 Gauss-Lobatto points in the pseudo-radial direction has no 
significantly effect on the Poincare surface of section. Finally, 
the Hamiltonian conservation is used as an independent accu- 
racy test on the computation. 

3. Acoustic ray dynamics in rotating stars 

In this section, we show that rotation strongly modifies the 
acoustic ray dynamics. Indeed, we find that, as rotation in- 
creases, the dynamics undergoes a transition from integrability 
(at zero rotation) to a mixed state where parts of the phase space 
display integrable behaviour while other parts are chaotic. 

The nature of a dynamical system is best investigated by con- 
sidering the structure of its phase space which contains both po- 
sition and momentum coordinates. We thus first introduce the 
Poincare Surface of Section (hereafter the PSS) which is a stan- 
dard tool to visualize the phase space (subsection 3.1). Then we 
describe the phase space structure at zero rotation (subsection 
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3.2) and the main features of the generic phase space structure at 
high rotation rates (subsection 3.3). The detail of the transition 
to chaos as rotation increases is analyzed in subsection 3.4. As 
this last subsection makes use of several specific tools and the- 
orems of dynamical systems theory, it might be skipped at first 
reading. 



3.1. Phase space visualization : Tine Poincare surface of 
section 

As shown in subsection 2.4, acoustic rays with a given L, are 
governed by a Hamiltonian with two degrees of freedom Hr. The 
associated phase space is therefore four-dimensional and diffi- 
cult to visualize. A PSS is constructed by computing the inter- 
section of the phase space trajectories with a chosen (2N - 1)- 
dimensional surface, where is the number of degrees of free- 
dom of the system. If H is time-independent, then energy con- 
servation implies that phase space trajectories stay on a (2N- 1)- 
dimensional surface. The PSS is thus a (2A^-2)-dimensional sur- 
face in general and a 2-dimensional surface in the present case. 




Fig. 1. (Color online) Intersection of an outgoing acoustic ray 
(red/arrow headed) with the r = fpiff) curve (magenta/thick). 
The point on the associated PSS is specified by the colatitude 
and the scaled latitudinal wavenumber component kg /lj at the 
intersection. 

Different choices are possible for the PSS although some 
conditions are required t o obtain a good description of the dy- 
namics (see for example lOtd I1993T) . First, in order to provide a 
complete view of phase space, the PSS must be intersected by 
all phase space trajectories. Here we chose the curve rp{6) - 
rs{0) - d situated at a fixed radial distance d from the stellar sur- 
face rs(0) displayed on Fig. [T] As shown in Appendix B.l for 
the non-rotating case, the distance d can be chosen such that all 
relevant trajectories intersect this curve. The second condition is 
that, given a point on the PSS, the next point on the PSS has to 
be uniquely determined. This is ensured by counting the inter- 
section with rp(ff) - rs(ff) - d only when the trajectory comes 



from on side of the r - rp(0) curve (here we consider the trajec- 
tories coming from the inner side). Finally, the coordinate sys- 
tem used to display the PSS is chosen in order that any surface 
of the PSS is conserved by the dynamics in the same way as 
four-dimensional volumes are preserved in phase space. The co- 
ordinates [6, kg] where kg is the latitudinal component of k in the 
natural basis (E^, E", E'^) associated with the coordinate system 
[( - rs{6) - r, 9, 0] fulfill this condition (as shown in Appendix 
B.2). 

PSS have been obtained by following many trajectories of 
different initial conditions. The number of trajectories together 
with the time during which they are computed determine the 
resolution by which the phase space is investigated. In princi- 
ple, we should display PSS computed for different values of the 
frequency cj. However, as u is varied in the range of frequency 
here considered, we found that the PSS remained practically un- 
changed. As discussed in subsection 2.4, this stems from the 
fact that the dynamics of the frequency-scaled wavevector kjtij 
is weakly dependent on w in this frequency range. 

3.2. The non-rotating case = 

The PSS at Q = is described in this subsection. It will serve 
as a reference to investigate the evolution of the dynamics with 
rotation. Due to spherical symmetry, the norm of the angular 
momentum with respect to the star center 



L = .\kl + 



sin0 



(20) 



is a conserved quantity. It follows that the intersection of any 
trajectory with the PSS belongs to a curve of the form: 



ICg - 



sinO 



(21) 



For L- = 0, these are the two straight lines kg = +L (see Fig. [3]) 
while Eq. ( 1211 1 yields a closed curve for L, the trajectories 
being constrained to latitudes smaller than arcsin(|L;|/L) (see 
Fig. |5]). This curve varies from a near rectangle to an ellipse as 
L. grows from 0^ to L. 

The simplicity of the PSS reflects the fact that the system 
is integrable (( |20l i indeed provides the second invariant (in ad- 
dition to H, ) of the reduced two-degree-of-freedom dynamical 
system). Every integrable system is structured in A^-dimensional 
surfaces which are associated with specific values of the con- 
stants of motion. This means that any trajectory is constrained to 
stay forever on one of these surfaces. They are called invariant 
tori because they are invariant through the dynamics and they 
have a torus-like topology. As we shall verify in the following, 
they play a crucial role in the transition to chaos as well as in 
the modes construction. The PSS at O = actually visualize the 
intersection of these tori with the r = rp surface. 

Importantly, the invariant tori can be of two different types 
which determine their fate once the rotation is increased. 
Rational (or resonant) tori are such that all trajectories on the 
torus are periodic orbits (that is trajectories that close on them- 
selves in phase space). In contrast, iiTational tori are such that 
any trajectory densely covers the whole torus. 

3.3. Phase space structure at high rotation rates 

The main features of the phase space at high rotation rates are 
shown in Fig. |2] where the PSS at Q = Q.59Q.K is displayed 
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Fig. 2. (Color online) PSS at Q = 0.59Q/S- and typical acoustic rays associated with the four main phase space structures : (a) 
a 2-period island ray (blue/dark grey) and the associated periodic orbit with endpoints a and b (orange/light grey), (b) a chaotic 
ray (red/grey), (c) a 6-period island ray (magenta/light grey) and (d) a whispering gallery ray (green/light grey). On the PSS, 
(colored/grey) symbols (diamonds for the chaotic and whispering gallery rays, crosses for the 2-period and 6-period island rays) 
specify the po ints where these trajectories cross the PSS. 



together with four acoustic rays shown on the position space 
as well as on the PSS. These rays belong to the three differ- 
ent types of phase space structures always present at high ro- 
tation rates. First, a large chaotic region appears (the grey region 
on Fig.|2]i. Chaotic rays, e.g. the red ray, are not constrained to 
stay on tori (that is on one-dimensional curves on the PSS) and 
thus fill up densely and ergodically a phase space volume (that 
is a surface on the PSS). Second, the island chain embedded in 
the large chaotic region is a common structure of phase space at 
high rotation rate. An important property of the island chain is 
to be structured by invariant tori centered on the periodic orbit 
of period 2 (the orange ray). The PSS also shows smaller island 
chains like the one formed around a 6-period periodic orbit (see 
the magenta ray). However, contrary to the 2-period island chain, 
such structures survive only up to a certain rotation rate. Third, 
the undulated curves present in the high kg region are formed 
by whispering gallery type trajectories (like the green ray), that 
is trajectories following the outer boundary. The associated tori 
correspond to the deformation of non-rotating tori which have 
not been destroyed yet at this rotation rate. Overall this type of 
phase space organization is typical of mixed systems with coex- 
istence of chaotic regions and invariant tori (the structures en- 
countered in integrable systems). 

We note also that the main phase space structures are dynam- 
ically independent since no trajectory can cross from one region 
to the other. We will show in section 4 that the very existence of 
these structures enables to organize the spectrum into indepen- 
dent frequency subsets. In the following subsection, the generic 
character of these structures is checked by computing the PSS at 
different rotation. 



3.4. Transition to cliaos 

The evolution of the dynamics with increased rotation is first 
described for L, = and then for 0. 



3.4.1. TheL- = case 

The PSS computed at the three rotation rates Q./Q.K = 
[0,0.15,0.32] respectively corresponding to the three flat- 
ness e - [0,0.01,0.05] are displayed in Fig. [3] to illus- 
trate the effect of a small centrifugal deformation on the ray 
dynamics. This perturbation of the integrable Q = sys- 
tem is very similar to one described by the KAM-theorem 
(Chirikov, 1979; Giannoni et al.', '1991'; 'Ott, 1993; GutzwiuiB, 
1 990; Lichtenberg & Lieberman. 1992;.Lazutkin, 1993 ). Indeed, 
for a smooth small perturbation of an integrable Hamiltonian, 
this theorem states that the tori of the integrable system which 
survived the perturbation occupy most of the phase space vol- 
ume. More specifically most of the irrational tori while being 
continuously perturbed can still be associated with N invariants, 
thus keeping their torus structure. This is the case of the undu- 
lated lines observed in the high kg domain of Fig. |3] By contrast, 
all rational tori are immediately destroyed for a non-vanishing 
perturbation. The KAM theorem implies that, despite the fact 
that the destroyed rational tori form a dense set in the phase 
space, the volume they occupy vanishes as the perturbation goes 
to zero. 

The theory of KAM-type transition to chaos also describes 
how resonant tori disappear. In our case, they correspond to 
one-dimensional curves on the PSS, formed by families of pe- 
riodic orbits. All orbits of the same torus will have the same pe- 
riod q. The so-called Poincare-Birkhoff theorem predicts that a 
(smooth) small perturbation will transform this one-dimensional 
curve into a chain of q stable points belonging to the same pe- 
riodic orbit and surrounded by stable islands, intertwined with 
q unstable periodic points. A small chaotic zone appears in the 
vicinity of the unstable fixed points and grows with the perturba- 
tion. The stable islands have themselves an intricate self-similar 
structure of small island chains surrounded by invariant structure 
(tori). This phenomenon is illustrated at Q./Q.K - 0.15 in Fig.|3] 
where, near the kg = curve, we can observe the 2-period is- 
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Fig. 3. Three L~ - PSS at small rotation rates showing the 
transition to chaos. The unit of koluj is uf^ . Island chains and 
chaotic regions appear around respectively stable and unstable 
periodic orbits. On the f2 - PSS, the straight lines correspond 
to intersections with mode-carrying-tori specified by the degree 
and radial order of the mode. 



land chain around the q-2 stable periodic points and the small 
chaotic region around the corresponding unstable points. This 
results from the destabilization of the resonant torus associated 
with the periodic orbits along the diameters of the spherical star 
Note that the width of the island chains (resonance width) is ex- 
pected to be approximately proportional to the square root of the 
perturbation, and decreases with q. 

What occurs for large perturbations following the KAM-type 
transition of integrable Hamiltonians has been studied in many 
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Fig. 4. Three L, = PSS visualizing the evolution of phase space 
as a function of rotation. All these PSS show the 2-period peri- 
odic orbit island chain embedded in a central chaotic region. As 
rotation increases, the 2-period island chain moves towards the 
equator while the central chaotic region enlarges. Note that the 
first Q = 0.32Q/i' PSS is displayed with a diff'erent scale in Fig. [3] 
The units are the same as in Fig. [3] 



systems. The general phenomenology which emerges also corre- 
sponds to what we observe in our system for increased rotation 
(see the PSS of Fig. g] computed for Q/Qj^ = [0.32,0.59,0.81] 
corresponding to the flatness e - [0.05, 0.15, 0.25]). The surviv- 
ing irrational tori as well as the island chains are progressively 
destroyed. This leads to the enlargement of the chaotic regions 
which were originally confined by these tori. This is illustrated 
in Figs. [3] and |4] where the surface of the central chaotic region 
becomes progressively larger with rotation. The island chains 
typically undergo a series of bifurcations for increasing pertur- 
bation; the most common bifurcation is the period-doubling one, 
where a stable periodic orbit of period q is changed to an unsta- 
ble orbit plus a stable orbit of period 2q. As the sequence of bi- 
furcations goes on, stable orbits are of longer and longer period 
until they eventually disappear. The destruction of stable regions 
is visible between Q - 0.59Qa: and Q - O.SlQjf (Fig.©, as the 
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Fig. 5. Three L, = OA/ljo PSS visualizing the evolution of phase 
space as a function of rotation. The evolution is qualitatively 
similar to the Z,, = case although the enlargement of the cen- 
tral chaotic region occurs at higher rotation rates. Same units as 
in Fig. [3] 



6-period island chain embedded in the chaotic central region at 
D. = 0.59Q.K has disappeared a higher rotation. As mentioned 
above, the fact that the largest stable island originates from a 
short periodic orbit (here a 2-period periodic orbit) is also a com- 
mon feature of the KAM-type transitions to chaos. 

While not visible on this figure, a zoom on other regions of 
the PSS would show the same process going on at small scales. 
It is however clear that the irrational tori associated with large 
values of L survive lo nger, this proper ty being also encountered 
in classical billiards dLazutkinl il99l . where tori close to the 
billiard boundary are the most robust. 



Qualitatively, the transition to chaos is very similar to the -0 
case. This is shown in Fig. |5] where PSS computed for L, = 
0.4/(x)o, are shown for increased rotation. The main effect of in- 
creasing Lr is to delay the transition towards chaos to higher 
rotation rates. Indeed by comparing PSS computed at the same 
rotation rate (see Fig.|4]and Fig.jS), one observes that the central 
chaotic region is more constrained by surviving tori for larger 

values. For example at f2 = 0.320.^ the central chaotic re- 
gion is much more developed for L, = than for L, = 0.4/aJo. 
The Q - Q.59Q.K PSS provides another example since for 
L, = 0.4/aJo the island chain associated with the 6-period orbit 
is separated by a surviving KAM tori from the central chaotic 
region while such a stable structure has already been destroyed 
for L. - 0. Finally, at Q = Q.SIQ.K, we can observe that the 
central chaotic region for L. = 0.4/(^0 contains visible surviving 
structures while this is not the case for L, - 0. 

The slower transition to chaos can be interpreted as be- 
ing due to the angular constraint - arcsin(|L-|/L) < 6 < 
arcsin(|L,|/L) imposed on the dynamics. This is compatible with 
the fact that for infinite L, the trajectories are confined to the 
equatorial plane and the dynamics becomes integrable because 
of the circular symmetry of this asymptotic limit. 

4. The asymptotic properties of tlie acoustic modes 

In this section, we show that ray dynamics provides a qualitative 
and quantitative understanding of the high-frequency acoustic 
modes. 

The question to be addressed is that of how to construct 
modes, that is standing waves, from the short-wavelength propa- 
gating waves described by ray dynamics. Such modes construc- 
tion is expected to be approximately valid in the asymptotic 
regime of high-frequencies (this asymptotic regime is called the 
semi-classical regime in a quantum physics context). As men- 
tioned in the introduction, the answer depends on the nature 
of the Hamiltonian system. For integrable systems, each phase 
space trajectory remains on an invariant torus and this enables 
to construct modes from a positively interfering superposition 
of these travelling waves on the torus. This is not any more the 
case for chaotic systems, where the ray dynamics provides no 
invariant structures on which to build modes. 

Thus for integrable systems, the modes and the frequen- 
cies can in principle be explicitly determined from the acous- 
tic rays through well-known formulas called Einstein-Brillouin- 
Keller (EBK) quantization after the name of it s main contrib- 
utors. We shall recall the results obtained by Gough ( 1993|) 
when applying the EBK method to spherical stars (sub-section 
4.1). While this procedure is not applicable to chaotic sys- 
tems, other concepts and methods have been developed and 
tested in quantum physics to interpret the non-integrable dy- 
namics. These concepts have been also applied to other wave 
phe nomena such as t h ose ob served in e.g. micr owave resonators 
(St ockmann & Steinl 1 19901). lasing cavities dNockel & Stonel 
T997I) . quartz bloc ks (Elle gaard et al.L 1 19961) and underwater 
waves dBrown et al. , 2003). Their potential interest for helio- 
seismology has been suggested, although no t demonstrated, by 
lPerdangl Tl988h and iKosovichev & Perdand (l9m . Here, we 
shall apply them to the non-integrable ray dynamics of rapidly 
rotating stars. More specifically, we have seen in section 3 that 
the ray dynamics of such stars corresponds to a mixed sys- 
tem where parts of phase space display integrable behaviour 
and other parts chaotic dynamics. In such case, the organiza- 
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tion/classification of modes in the semi-classical regime is ex- 
pected to closely follow the structure of phase space. Near- 
integrable regions of phase space like the island chains are 
amenable to EBK quantization, leading to a regular frequency 
spectrum, while the modes originating from chaotic regions have 
an irregular frequency spectrum with generic statistical proper- 
ties. Another important information provided by ray dynamics 
is the averaged number of modes that can be constructed from a 
given phase space region, this number being proportional to the 
volume of the region considered. 

In the following, we explicit these concepts and methods 
in the context of stellar acoustics (sub-sections 4.1, 4.2, 4.3). 
Then, their relevance to describe the asymptotic properties of 
the acoustic modes is tested by confronting their predictions 
to the actual properties of high-frequency acoustic modes (sub- 
sections 4.4, 4.5, 4.6, 4.7). These modes are axisymmetric modes 
in the frequency range [9tiJi, 12<yi], being the lowest acoustic 
mode frequency of the stellar model considered. They have been 
computed for a Q = 0.59Qa: uniformly rotating = 3 polytropic 
model of star and under the same assumptions as for ray dynam- 
ics (adiabatic perturbations, no Coriolis acceleration. Cowling 
approximation). 

4.1. The integrable case Q.-Q 

To build modes from the ray dynamics, the wave-like solution 
*P = A(x) exp[/0(x)] is regarded as a function of the phase space 
variables x,k - that is subsequently projected onto the posi- 
tion space. The condition that exp[/(l)(x)] be single-valued on the 
position space requires that for any phase space trajectory that 
closes on itself in the position space, the variation of O along 
this closed contour is a multiple of In. As trajectories of an in- 
tegrable system stay on invariant tori, this condition leads to the 
EBK formula: 

r k-dx ^2n{ni+—) (22) 
Jcj 4 

where Cj is any closed contour on a given torus and rij and y6y 
are non-negative integers. The integer called the Maslov in- 
dex is introduced to account for a 7r/2 phase lag that must added 
each time the contour crosses a caustic. Indeed, the caustic cor- 
responds to the boundary of the torus projection onto position 
space; the amplitude A taken in the position space is discontinu- 
ous t here, leading to the 7r/2 phase loss (see lKeller & Rubinowl 
[mOii for details). 

While this expression is valid for any closed contour on the 
torus, it can be shown that it gives the same condition for all 
contours that can be continuously deformed to the same one. 
Thus in fact the EBK quantization yields independent con- 
ditions, as only topologically independent closed paths exist 
on a A^-dimensional torus. As these closed paths do not need to 
be actual trajectories of the dynamical system, the usual way to 
construct EBK solutions is to choose contours for which the for- 
mulas are simple to compute. The quantization conditions thus 
select a particular torus on which a mode can be built, irrespec- 
tive of the fact that the torus is resonant or non-resonant. For 
spherical stars, three independent contours on a torus specified 
by L, and w can be obtained by varying one of the spheri- 
cal coordinate s and fixing the other two. Using similar contours, 
lGoughl(II993h obtained the three conditions: 

J [ ^2 - 72 ) " - 2^''' -^ = ^+5' ^^^^ (23) 



where n, £, m are non negative integer, r, and re being the in- 
ternal and external caustic respectively. Note that L - a>L and 
L, = cjL^. The associated eigenmodes have been also explic- 
itly constr ucted from the trajectories on the selected torus. As 
shown by iGoughl (1 1 9931) . the result of this EBK quantization 
is practically identical to the usual asymptotic theory of acous- 
tic modes in spherical stars which uses the separability of the 
three-dimensional eigenvalue problem and a WKB approxima- 
tion of the resulting 1-D boundary value problem in the radial 
direction (in the usual analysis, L = V^(7TT), differs from the 
EBK result, especially at low £ values). This shows that the in- 
tegers n, {, m derived from the EBK quantization conditions do 
correspond to the order, degree and the azimuthal number of the 
acoustic modes in spherical stars. 

The tori on which the eigenmodes are constructed can be vi- 
sualized on the PSS. For example, the (f, n, m) = (8, 10, 0) mode 
is associated with the torus L - +{£ + 1 /2)/w„/_„„ L. - which 
imprint on the PSS are the straight lines kg - ±L. The inter- 
section of various mode-carrying tori with the L, - Q PSS are 
shown in Fig. |3] High radial order modes approach the central 
^fl = line while large degree modes occupy the high kg region. 

4.2. Chaotic systems 

It has been widely recognized in the past decades that most dy- 
namical systems are not integrable and therefore display various 
degrees of chaos. The quantum chaos field has studied quan- 
tum systems whose short-wavelength classical limit displays 
such c haotic behavior As has been recognized early by Einstei^ 
(Il917l) the EBK quantization explained in the above paragraph 
cannot be applied to such systems. Indeed, no A^-dimensional 
invariant structure exists on which to apply conditions of con- 
structive interference like the EBK condition. Rather, the semi- 
classical limit of these chaotic systems yields a Fourier-like for- 
mula which connects the set of all classical periodic orbits to 
the w hole spe ctrum. This formula, called Gutzwiller trace for- 
mula dOutz wi lier, 1990) is much more delicate to use than EBK 
formulas, since it represents a divergent sum from which infor- 
mation can be extracted only through refined mathematical and 
numerical methods. 

On the other hand, the very complexity of chaotic systems 
leads to statistical universalities. Indeed, in a similar way as sta- 
tistical physics emerges from the random behavior of individual 
particles, in chaotic systems the randomness induced by chaos 
leads to robust statistical properties of eigenmodes. Contrary to 
modes of integrable systems, which are localized on individual 
tori selected by the EBK conditions, in chaotic systems modes 
are generally not associated to a specific structure in phase space, 
and are ergodic on the energy surface. It has been found that 
one can model such systems by replacing the Hamiltonian by 
a matrix whose entries are random variables with Gaussian dis- 
tributions. Such ensembles of Random Matrices, which contain 
no free parameter but take into account the symmetries of the 
system, can give precise predictions, which have been found 
to describe accurately many statistical properties of the modes 
of systems with a chaotic classical limit. This has b een conjec- 
tured and checked numerica lly for many systems CBohigas et al.L 
il984tlGiannoni et alill99ll) . 

The comparison with Random Matrix predictions is often 
done through the statistical analysis of the frequency spectrum. 
In general the density of modes as a function of the frequency 

d(cj) = J]6iaj-cj„) (24) 
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where 6 is the delta function, can be written as the sum of two 
functions 



(25) 



The quantity d"^'(a)) (hereafter called the Weyl term) is a smooth 
function which describes the average density of modes at a given 
frequency. It has been known from the beginning of the twentieth 
century that this term can be calculated from general geometric 
features of the system such as phase space volume and therefore 
is independent of the chaotic or integrable nature of the dynam- 
ics (d'"\cL>) is estimated for stellar acoustic modes in subsection 
4.7). In contrast, the function d^'"" which describes the fluctua- 
tions around the mean position of eigenfrequencies strongly de- 
pends on the nature of the dynamics. (Note that most textbooks 
on this subject use the quantum physics terminology, that is 'en- 
ergy level' instead of frequency and 'density of states' instead of 
density of modes). 

The spectra of integrable systems are predicted to be un- 
corrected, and in gener al this leads to fluctu ations given by 
the Poisson distribution (lBerrv&Taboiill977h . In contrast, for 
chaotic systems these fluctuations should be given by Random 
Matrix Theory. The theory of Random Matrices has therefore 
been developed in order to compute analytically the predic- 
tions for specific quantities, which in turn could be compared 
to numerical data for real systems. A popular quantity to de- 
scribe fluctuations in spectra is the spacing distribution Pi6), 
which is the distribution of spacings in frequency between con- 
secutive eigenfrequencies, once the frequency differences have 
been rescaled by the Weyl term such that the average spacing 
is one. The function P(6) measures the coiTelations at short dis- 
tances in frequency in the spectrum, and therefore does not give 
information about all statistical properties, but is nevertheless 
very useful since the predictions are strikingly different for the 
Poisson distribution and for Random Matrix Theory. While the 
Poisson distribution coiTesponds to P{6) = exp(-(J), the pre- 
diction of Random Matrix Theory is the Wigner distribution 
P{6) = 7t6/2 exp{-7T6^ /4), which displays frequency repulsion 
(level repulsion in the quantum terminology) at short distances 
(small 6) and falls off faster than Poisson at large 6. 

4.3. Mixed systems 

We have seen in section 3 that the acoustic ray dynamics in ro- 
tating stars has a mixed character as chaotic regions coexist with 
stables structures like island chains or invariant tori. Such mixed 
system are actually the most common in nature, completely in- 
tegrable and chaotic systems being limiting cases. 

In the context of quan tum chaos, seminal studies of these 
systems bv lPerciva]|(ll973h and lBerrv & RobnikI ( 11984 enabled 
to conjecture that a good description of their spectrum at high 
energy is obtained by quantizing independently the structured 
and chaotic parts. While a zoom on island regions would reveal 
a complex structure involving chaotic trajectories and chains of 
smaller islands, these small scale details can be neglected for 
the island quantization if the mode wavelength remains large 
as compared to these scales. Instead, the presence of a large 
number of invariant structures constrains enough the dynamics 
to make the system similar to a purely integrable structure to 
which EBK quantization applies. These region are then called 
near-integrable. 

Thus subsets of modes can be associated to the different 
near-integrable island chains, while other subsets correspond to 
the chaotic zones. In each subset, the modes behave as if they 



were constructed from an isolated system. Thus in mixed sys- 
tems the frequency spectrum can be described as a superposition 
of independent frequency subsets associated with the different 
phase space regions. Subsequent works have shown this picture 
to be a good approximation of actual spectra, although in some 
cases certain correlations are present between the frequency sub- 
sets due to the presence of partial barriers in phase space or to 
the existence of mod es localized at the border between zones 
(iBohigas etal.Lll993h . 

The acoustic ray dynamics of rapidly rotating stars being of 
this mixed type, one can expect such an organization of the spec- 
trum to be valid, even though the modes are of quite large wave- 
length compared to previous studies in quantum chaos. To test 
this hypothesis, it is convenient to compute a phase-space rep- 
resentation of the modes. Indeed, the chaotic or near-integrable 
zones are well-defined in phase space while their projections in 
position space, were modes are usually pictured, are generally 
much more difficult to separate. 

4.4. Associating modes to rays 

Constructing phase space representations for modes has been 
envisioned since the beginning of quantum mechanics, since 
it enables to test the quantum-classical correspondence accu- 
rately. Contrary to states of a classical system which are de- 
fined by a point in phase space, modes have always a finite ex- 
tension in phase space since they have a finite wavelength and 
their localization in wavenumber space is, according to Fourier 
analysis, inversely proportional to their localization in physi- 
cal space. Any mode occupies a finite volume of the order of 
(In)^ in the physical/wavenumber phase space (a (In/cj)'^ vol- 
ume in the physical/scaled-wavenumber (x, k) phase space or a 
(IttH)^ volum e in the p osition/momentum phase space of quan- 
tum physics). IWignerl Il932) was the first to construct a phase 
space function representing a mode, but this so-called Wigner 
function has the disadvantage of being positive or negative, and 
therefore cannot be interpreted as a probability distribution of 
the mode. To circumvent this problem, one way is to smooth the 
Wigner function by a Gaussian convolut ion. The resulting fu nc- 
tion, called Husimi distribution function dChang & ShiLll986l) . is 
always real and nonnegative and can be understood equivalently 
as the projection of the mode onto a Gaussian wave packet cen- 
tered around x and k: 



Wx, k) = 



"Vix') exp(-||A:' - Jc||^/(2A^)) exp(/wjt ■ x')dx' 



(26) 



where ^'{x') is the mode and exp(-||j[:' - x\\^ / (2A^^)) expiicok) 
the Gaussian wavepacket. In this expression, the width of the 
wavepacket A^^ determines the resolution of the Husimi function 
in the spatial direction, the resolution in the scaled-wavenumber 
Aj being such that A^^A^; x l/oj. These quantities determine the 
minimal extent of the mode representation in both directions. 

The computed modes are 3D modes and they shall be com- 
pared with the reduced ray dynamics computed on a 2D merid- 
ional plane. As shown in the spherical case by Gough ( 1993), the 
amplitude of a 3D axisymmetric mode constructed from acous- 
tic rays obtained on neighboring meridional planes decreases 
as (r sin0)*~'^^' because away from the rotation axis the dis- 
tance between the planes and thus the density of rays increases. 
Thus, the computed 3D modes have been scaled by (rsinfl)''^^' 
to better represent the mode amplitude on a meridional plane. 
Moreover to obtain a phase-space representation limited to the 
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Fig. 6. (Color online) Four axisymmetric modes and their phase 
space representation : (a) a 2-period island mode (blue/dark 
grey), (b) a chaotic mode (red/grey), (c) a 6-period island 
mode (magenta/light grey) and (d) a whispering gallery mode 
(green/light grey). The amplitude distribution of the scaled mode 
is represented by its nodal lines (black) and its positive 
(full colored/grey lines) and negative (dashed colored/grey lines) 
level curves. The Husimi distributions of the four modes com- 
puted for the same Aj = OAIR^ are represented by their level 
curves. 

PSS, we actually computed the Husimi's distribution function of 
the ID cut of the mode taken along the PSS: 



"H(sJ) = 



r(s')exp(-(s' - sf/(2A^^))expiiojh')ds' 



(27) 



where *?' is the scaled version of the mode - Pja solution of 
Eqs. (|6l), (|7]i and ([8]), 5 is a curvilinear coordinate along the curve 
r - Tp and k - k ■ Cp is the scaled-wavenumber in the direction 
tangent to this curve, Cp being a unit vector tangent to the curve. 

is then expressed in terms of the PSS coordinates using the 
following relations between [0, kg] and [s, k]: 



1/2 

d0' k^kgE^-e^ 



(28) 



The vector is defined in appendix B.2. The integral dZTl l is 
performed in the interval [0-jT,9+n], the mode being prolonged 
by symmetry outside the [0, n] interval. 

The Husimi function has been computed for the axi- 
symmetric modes of the Q. = 0.59Q.K star and its contour plot 
is compared with the PSS of the ray dynamics computed for the 
same star model. Fig. |6] illustrates this process by showing the 



position space as well as the phase space representation of four 
typical modes. As can be observed, the modes can be clearly 
associated with one of the main structure of the phase space 
namely, the 2-period island chain, the large central chaotic re- 
gion, the 6-period island chain or the whispering gallery region. 
On the PSS, we note however that the Husimi function is sym- 
metric with respect to kg while the dynamics is not. This dif- 
ference is due to the fact that the PSS is constructed only with 
kr > intersecting trajectories while the Husimi function com- 
puted from the mode cut on the r - rp contains no information 
about the sign of kr- Nevertheless, in the high-frequency inter- 
val {9(jL>\, \2(jj\\ that we studied in detail, any ambiguity on the 
phase space location can always be resolved using the additional 
information on the mode distribution in the position space. In 
this frequency interval, we thus classified the modes according 
to their localization in phase space distinguishing the 2-period 
island modes, the chaotic modes, the 6-period island modes and 
the whispering gallery modes associated with the correspond- 
ing phase space regions. As a result, the full frequency spectrum 
can be discomposed into sub-spectra associated with phase space 
structures. Fig.|7]displays the four sub-spectra. 

In the following, we shall analyze these sub-spectra and test 
whether the Percival and Berry-Robnik conjecture described in 
subsection 4.3 applies to acoustic modes in rapidly rotating stars. 
We first study the regular character of the sub-spectra issued 
from near-integrable phase space regions (sub-section 4.5) and 
then consider the spectrum of chaotic modes (sub-section 4.6). 

4.5. The regular spectra 

A spectrum is said to be regular if it can be described by a func- 
tion of integers, being the degree o f freedom of the system . 
In accordance with pr evious stud ies by [Li gnieres et al.l (l2006l) . 
iLignieres & Georgeot! fcoOS) and lReese et al. (.2008;), the spec- 
trum of the 2-period island modes is well-fitted by the empirical 
formula 



a)„c - n6„ + 16 c + a 



(29) 



confirming the regular nature of this spectrum that is also clearly 
apparent on Fig. [Tja). 

The 6-period island mode spectrum shown in Fig.|7tc) is also 
regular as it is closely fitted by the even simpler formula: 



n 6' + a' 



(30) 



Indeed, the root mean square error between this empirical fit and 
the actual spectrum is equal to 1 .9 percent of 6]^ (where 5', has 
been determined as the mean of the spacing between consecu- 
tive frequencies and a' is fixed such that the model is exact at a 
reference frequency). 

While a simple linear law such as Eqs. ( |29l ) or ( [30l ) does 
not apply to the whispering gallery modes, there are strong evi- 
dences that this sub-spectrum is also regular. Thanks to the reg- 
ularity of the nodal lines pattern (as apparent in Fig. |6jd)), two 
integers corresponding to the number of nodes along the polar 
axis and to the number of nodes following the internal caustic 
can be easily attributed to each modes. When plotted against the 
number of caustic nodes (as in Fig. |2ld)), the spectrum clearly 
shows a regular behavior. The fact that the function of these 
two integers describing the spectrum is however not as simple 
as Eqs. ( |29] l or dSOl l is expected from what we know about the 
regular s pectrum of high degree modes in spherical stars (see for 
example lChristensen-DalsgaardI (Il980l) ). 
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The regularity of the three sub-spectra issued from near- 
integrable phase space regions is fully in accordance with the 
Percival's conjecture. An important consequence is that theoret- 
ical model of these spectra can in principle be obtained from the 
EBK quantization of the invariant structures of the correspond- 
ing near-integrable regions. As a result we should be able to re- 
late the potentially observable quantities (5„, 6( or to the star 
properties. In practice, the standard method is first to construct 
the normal forms around the central periodic orbit in order to 
describe the dynamics in the island, and then use the EBK quan- 
tization scheme to fi n d the asy mptotic formula for the modes 
jBohigas et all 119931: iLazutkinl [1993 ). While such a program 
is outside the scope of the present paper, we men tion below 
the result obtained in iLignieres & GeorgeotI (|2008[) for the 2- 
period island modes using an equivalent procedure, which may 
be more physically transparent, and extend it to the 6-period is- 
land modes. 

As already noted, the propagation of acoustic waves in our 
system is similar to the propagation of light in an inhomoge- 
neous medium, the role of the medium index being played by 
l/cj. The construction of standing wave solutions between two 
bounding surfaces has been investigated i ntensively in the con- 
text of the study of laser modes in cavities (Ko gelnik & Li[T966l) 
and it consists in applying the paraxial approximation in the 
vicinity of the optical axis. While generally applied to homo- 
geneous media, this app roximation can be extende d to the non- 
homogeneous case as in lPermitin & Smirnovl (11996) . Applying 
this fo rmalism to the acousti c mode s associated with periodic 
orbits, ILignieres & GeorgeotI ( l2008l) found a model spectrum 
equivalent to Eq. ( |29] l with 



5„ - 



and 6c -2 



^ Csdcr/w^ 
fa ^^^/^^ 



(31) 



where cr is the curvilinear coordinate along the periodic orbit 
and the integral is computed between the end points of the or- 
bit (these points are shown in Fig. |2| for the 2-period and the 
6-period periodic orbits being denoted (a,b) and (a',b') respec- 
tively). The quantity vv'(cr) in the expression of 6( describes the 
spreading of the wave beam in the direction transverse to the 
periodic orbit and verifies a differential equation which depends 
on the sound speed and its transverse derivative taken along the 
periodic orbit. Moreover the integers n and f correspond respec- 
tively to the number o f nodes in the directions paral lel and trans- 
verse to the orbit (seel Lignieres & GeorgeotI (l2008l) for details). 

When applied to the 2-period periodic orbit, this theoretical 
expression of 6„ yields a value 0.5635wo which differs from its 
empirical value by only 2.2 percent. 

The 6-period island mode spectrum can be modeled in the 
same way. In the frequency interval considered, these modes 
have a similar structure in the direction transverse to the central 
orbit and should therefore be associated with the same f value. 
The model spectrum has thus the same form as Eq ( [30l ) where 



S'n 



(32) 



As for the 2-period modes, we find that this theoretical value of 
6'n differs by only a few percents from the empirical determina- 
tion of 6'n - 0.186(Wo- 

These two examples show that ray dynamics can provide a 
quite accurate model of the near-integrable spectra in the rel- 
atively low frequency regime considered. Moreover the model 
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Fig. 7. Frequency sub-spectra of four classes of axisymmetric 
modes : (a) the 2-period island modes, (b) the chaotic modes 
antisymmetric with respect to the equator, (c) the 6-period is- 
land modes and (d) some whispering gallery modes. For the sub- 
spectra (a) and (d), the height of the vertical bar specifies one of 
the two quantum numbers characterizing the mode. 



of the 2-period island m ode spectrum remains reasonably 

accurate at lower frequencies ( ILignieres et all 20061) and can be 
extended to non-axisymmetric modes dReese et aUl2008l) . 



4.6. The chaotic modes 

A large subset of the frequencies correspond to modes localized 
in the chaotic zone of phase space (the chaotic modes). As we 
have seen in subsection 4.2, one should not expect regular pat- 
terns for this part of the spectrum. Rather, the chaotic charac- 
ter of the phase space should be reflected in specific statistical 
properties of the sub-spectrum, which should follow predictions 
from Random Matrix Theory. To test this predictions, we have 
studied the distribution of the consecutive frequency spacings 
6i = - (L)i of the chaotic modes. Fig.|8]shows the integrated 

distribution A^(A) = P{8)d5 (with spacings normalized by 
the mean level spacing within the chaotic subset, as should be 
done). The distribution is constructed from the two independent 
distributions obtained for the equatorially symmetric and anti- 
symmetric modes, corresponding to around 187 modes in total. 
Although the difficulty to solve Eqs. (|6l), (|7]i, ([8]) prevents us to 
reach as large frequency samples as can be obtained for other 
systems ( Bohigas et aU 1984). the numerically computed A^(A) 
agrees well with the Random Matrix predictions, and is clearly 
different from the Poisson distribution typical of integrable sys- 
tems. This result indicates that these modes, selected by means 
of the comparison between their localization in phase space and 
the ray dynamics, have indeed the frequency statistics expected 
for chaotic modes. We therefore think that this confirms the va- 
lidity of the ray model, and gives a strong evidence that wave 
chaos actually occurs in rapidly rotating stars. 

We note that the modes we identify as chaotic are located in 
the chaotic zone but cover only part of it, as for the example of 
Fig.|6l We think this is partly due to the relatively low frequency 
considered, which prevents ergodicity of the modes to be clearly 
visible. In addition, it is known that certain low-energy eigen- 
functions of chaotic systems called scars are conc entrated along 
short periodic orbits of the system dHelleii Il984t) . In this case. 
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Fig. 8. Integrated spacing distribution A^(A) of the chaotic 
modes (full line). The dashed line shows A^(A) for the Gaussian 
Orthogonal Ensemble of Random Matrix Theory while the dot- 
ted line shows the integrated Poisson distribution typical of inte- 
grable systems. 

rather than being ergodic, some individual modes are effectively 
localized in the vicinity of such orbits. This effect can create 
some subsequences of low-energy modes with regular frequency 
patterns, even if the distribution of frequency spacings follows 
the predictions of Random Matrix Theory. Precise investigation 
of this phenomenon in the context of stellar acoustic modes may 
be important. 

4. 7. Predicting the number of modes in each sub-spectrum : 
The Weyl formula 

In this subsection, we show that ray dynamics enables to esti- 
mate the number of chaotic and island modes present within 
a given frequency interval. This information is complementary 
to the regular/iiTegular properties of the associated sub-spectra 
shown in the previous section and it is crucial to build an asymp- 
totic model of the frequency spectrum. 

We have seen in subsection 4.2 that the density of modes 
as a function of the frequency d(oS) can be written as the 
sum of a smooth part d'"'{cL>) and an oscillatory part d^'""{(jS) 
(see Eq. (IZSTi). Weyl, at the beginning of the twentieth cen- 
tury, derived analytically an asymptotic expansion of d"^'{cj) 
( IWevll[T9T2l) . The leading term of the Weyl formula can be ob- 
tained from general principles. We have seen (subsection 4.4) 
that in average a mode occupies a (2n)^ volume in the physi- 
cal/wavenumber phase space. The averaged number of modes in 
a given phase space volume can thus be estimated as the volume 
of phase space available divided by (Ijif , the volume occupied 
by one mode. In the following, we shall first verify that the lead- 
ing term of the Weyl formula yields a reasonable estimate of 
the number of modes in the case of spherically symmetric stars. 
Then, we shall calculate the phase space volume of the chaotic 
and island regions and confront the result with those obtained 
from the numerical computation of modes at Q = Q.59D.K- 

The Hamiltonian formulation H' — u - -y/cfX^TTjJ is best 
suited for this calculation. The averaged number of modes below 
a given frequency w reads: 

(In)" 

where 'V(w) is the phase space volume coiTesponding to ener- 
gies smaller than w defined as 

^{cj}^ r dk^dx^ (34) 
where J[(io) is the phase space region defined as H'(x, k) < oj. 



For spherical stars, the dynamics can be reduced to the 
one-degree-of-freedom dynamics characterized by the reduced 
Hamiltonian - ^Jcsir)^ikj + L-/r^) + w^. Applying the 

above formula to H^'''\ the double integral 'V{a)) can be inte- 
grated over kf to give: 

which is thus the estimated number of modes of given L and 
L, with a frequency smaller than co. This estimation can be com- 
pared to the results of the usual asymptotic theory (see Eq. (|23]l). 
Accordingly, N'^P^iu) - n-l/2, where n is the radial order asso- 
ciated to the frequency a> but is also the number of modes below 
the frequency oj for fixed values of L and L,. We therefore con- 
clude that in the ID-spherical case the first term of Weyl's for- 
mula yields a reasonable approximation of the averaged mode 
density. 

In rotating stars, to estimate N{a>) the number of modes 
below o) for a given L,, we use the two-degree-of-freedom 

Hamiltonian h[. - a> - ^c^^ikj^ + L^/(rsin0)^) + oj^, and in- 
tegrate the 4-dimensional volume integral in the wavenumber 
directions to obtain: 

47rijs,„ cj (r sin 0)2 

where S„, is the portion of the meridional plane surface where 
the integrand is positive. While providing the total number of 
modes, this expression does not give the fraction of chaotic, is- 
land or whispering gallery modes which are important quantities 
to model the spectrum. These quantities can nevertheless be ob- 
tained by computing the corresponding phase space volumes and 
by applying Eq. (33[ . This has been done for Q./Q.K = 0.59 since 
at this rotation rate we can compare the results of Weyl's formula 
to the numbers of chaotic and island modes obtained from the di- 
rect mode computation and the classification through the Husimi 
distribution. 

The 4-dimensional phase space volumes have been evaluated 
using a Monte-Carlo quadrature method: points are randomly 
chosen in a known volume Vm that includes the volume V to be 
computed. The proportion of points inside V approximates the 
ratio V/Vm, thus providing an approximate value of V. To de- 
cide whether a given point is inside or outside V, we used space 
filling trajectories on the torus delimiting the volume V. Two 
phase space volumes have been computed at O/Qj^ = 0.59. The 
first one includes the large chaotic region as well as the island 
chains around the 2-period and 6-period orbits. The second vol- 
ume corresponds to the 2-period island chain. The details of the 
calculation and an estimation of the error on the volume deter- 
mination are given in Appendix C. 

As a result, the leading term of the Weyl formula yields 
34 ± 2 modes in the 2-period island chain and 270 + 8 modes 
outside the whispering gallery region in the frequency interval 
[9.42wi, 11.85wi]. This value has to be compared with the 50 
island modes and the 276 modes outside the whispering gallery 
region obtained using the Husimi phase space representation of 
the modes computed in the same frequency interval. 

The difference between the estimation given by ( [33] ) and the 
actual number of modes in each subset of the frequency spec- 
trum corresponds most likely to the next order in the asymptotic 
expansion of the density of modes. Indeed, Eq. ( [33T l is only the 
first term in an asymptotic expansion, the next term being usu- 
ally proportional to the length of the boundary between phase 
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space zones. At relatively low frequency, this term can be signif- 
icant. Another source of imprecision can stem from the fact that 
at such relatively low frequency some partial barriers in phase 
space can trap island-like modes in the vicinity of the island, en- 
larging the effective size of the regular zone. Indeed, for some 
of the modes classified as island modes, the Husimi distribution 
is not entirely inside the island, the outer part remaining close 
to the island. Nevertheless, our study show that Eq. ( l33T l gives 
a reliable estimate of the relative sizes of the frequency subsets 
which can be obtained without any knowledge of the spectrum 
itself. 



5. Discussion and applications to asteroseismology 

In this section we discuss the validity of the assumptions of our 
asymptotic analysis, and the implications of our results for the 
asteroseismology of rapidly rotating stars. 

5.1. Assumptions of the asymptotic analysis 

The WKB assumption underlying the asymptotic analysis is not 
justified for the lowest frequency acoustic modes. While deter- 
mining the limit of validity of the different results presented 
here is outside the scope of the paper, we know that the reg- 
ularities of the 2-period island mode sub-spectrum are rele- 

vant down to about th e 5*^ ^ radial order acoustic pulsations (see 
iLipnieres etsl\ (HOO^ and' Reese etsl\ (|2008) for details). 

Another important precaution when applying the present 
analysis stems from the fact that the interpretation of the ray 
dynamics depends on the frequency range considered. Indeed, if 
extremely high frequency modes were to exist, their properties 
would closely follow the phase space structure up to its smallest 
details. This means for example that such modes could be as- 
sociated with the very small chaotic regions which exist inside 
the island chains or in between the surviving KAM tori of the 
whispering gallery regions. On the other hand, finite-wavelength 
effects have to be taken into account when interpreting the ray 
dynamics at finite frequencies. For example, the regularities of 
the island mode spectra in the [9(l)] , I2a)[] interval show that the 
small chaotic zones within the island chain can be overlooked 
in this frequency range. The same reasoning holds if one wants 
to interpret the ray dynamics at small rotation rates (see Fig. [3]). 
The tiny chaotic regions predicted by the KAM theorem could 
be interpreted as a proof for the existence of a chaotic mode 
frequency subset at vanishingly small rotation. However, these 
modes should have such a short wavelength to "fit in" the chaotic 
region that they may simply not exist (because they are strongly 
dissipated by diffusive effects or their frequency is so high that 
they are not reflected at the surface). 

Ray dynamics can not account directly for the coupling ef- 
fects between two modes associated with two dynamically iso- 
lated regions of phase space (as occurs for example in the well 
known tunneling effect). Indeed, while trajectories can not cross 
the dynamical barrier between the chaotic and the island chain 
regions, an island mode can be present on either side of the bar- 
rier if its frequency is very close to the one of a chaotic mode 
(and vice versa). As usual for mode avoided crossings, the mode 
distribution can thus be significantly perturbed by the coupling 
but the frequency is only slightly affected especially in the high- 
frequency regime. Quantifying the effects of such avoided cross- 
ings would require a specific study. 

The WKB assumption that the wavelength is much shorter 
than the typical background lengthscale breaks down for low- 



frequency acoustic mode but also when the typical lengthscale 
of the stellar model becomes very small. This can occur in real 
stars, especially at the upper limit of the core convective zone 
where strong composition gradients built in during evolution. 
The effect of such a discon t inuity has been studied for spher- 
ical stars jVorontsovl 119881: iGoughl 1 19901: iBallot et all |2004|) 
and has been found to add an oscillatory component to Tassoul's 
asymptotic formula, but not to remove the asymptotic structure 
altogether. To treat properly the discontinuities in non-spherical 
stars, the ray dynamics approach has to be extended by taking 
into account the splitting of rays at the discontinuity, coiTespond- 
ing to the reflected and transmitted waves. Once this is incorpo- 
rated in the formalism, quantum chaos techniques can be applied 
as was done for billiards in Bluemel et al. ( 1996ah and the Weyl 
formula can also be computed tPrange et al.i . il99 o). 

Apart from the treatment of eventual sound-speed discon- 
tinuities, the asymptotic analysis presented in this paper for 
uniformly rotating polytropic stellar models can be readily ap- 
plied to realistic stellar models. The details of the dynamics will 
change as they depend on the sound speed distribution of the 
model considered. However, we do not expect that the mixed 
character of the dynamics and thus the irregular/regular nature 
of the spectrum will change. This has to be confirmed by specific 
ray dynamics studies. In particular, the effects of the advection 
by a differentially rotating flow should be investigated. 

Lifting the two assumptions concerning the adiabaticity of 
the perturbations and the Coriolis force omission should not sig- 
nificantly modify the results of the asymptotic analysis. Non- 
adiabatic calculations are known to have a small effect on the 
frequency values while the legitimate omission of the Coriolis 
force for high frequency motions is a l ready relevant for quite 
low frequency as shown in lReese et al. I (l2006h . 

5.2. Implications for mode identification 

The asymptotic analysis provides qualitative and quantitative in- 
formations which can be used to identify high frequency acous- 
tic modes in an observed spectrum. First the basic structure of 
the spectrum can be readily deduced from the ray dynamics 
phase space structure visualized by the PSS. Indeed, we have 
seen that the Q. - 0.59Qjf PSS correctly predicts that the spec- 
trum of axi-symmetric modes is a superposition of four fre- 
quency subsets, three regular and one irregular. If we now look at 
the Q - PSS, we see that the spectrum structure should 

be similar except that the regular sub-spectrum associated with 
the 6-period island chain is no longer present since this island 
chain has disappeared at this rotation rate. In the same way, at 
Q. = 0.15f2/i-, we expect a simple superposition of a whispering 
gallery sub-spectrum and a 2-period island mode sub-spectrum 
since the chaotic regions are not sufficiently developed at this ro- 
tation rate. Such informations, while only qualitative, are crucial 
to guide the identification process. Moreover this information is 
obtained at a relatively low computing cost since the PSS calcu- 
lation is much-less demanding than the numerical computation 
of modes and frequencies. 

Then, the EBK quantization of the near-integrable regions 
provides the values of the uniform frequency spacings (as given 
by Eq. ( |29l )) that should be present in the observed spectrum. 
When analyzing an observed spectrum, the star model is not 
known, thus only estimates of these uniform spacings can be ob- 
tained. However, these estimates enable to focus the search for 
regularities on a smaller range of values. 

Finally as we also know the frequency statistics of chaotic 
modes and the number of modes in each sub-spectrum (through 



F. Lignieres and B. Georgeot: Asymptotic analysis of high-frequency acoustic modes in rapidly rotating stars 



15 



Weyl's formula) the asymptotic analysis enables to construct 
asymptotic spectra. The chaotic mode frequencies can be ob- 
tained as a realization of the Wigner distribution although, in 
this case, their frequencies could not be individually identified 
with observed ones. Nevertheless, such a synthetic spectrum 
should be very useful to test identification methods, especially 
the search for the regularities hidden in the complete spectrum 
(see below). 

Among the additional informations which can help con- 
straint the mode identification are the mode visibility and the 
mode excitation. The excitation mechanism have been studied 
so far in the spherically symmetric case and need to be recon- 
sidered for rapidly rotating stars. The mode visibility also de- 
serves a specific study notably the calculation of the intensity 
variations induced by the oscillation in a gravity darkened at- 
mosphere. However, the visibility strongly depends on the can- 
cellation effects on the disk-integrated light. Here, we can esti- 
mate this effect by integrating the surface Lagrangian tempera- 
ture perturbation of the axisymmetric modes computed for the 
£1 - 0.59Qa: rotating polytropic star The disk-averaging factor 
is defined as: 



D(i) 



6T{0, 0)dS ■ Ci 



(37) 



where / is the inclination angle between the line-of-sight and the 
rotation axis, Ci is a unit vector in the observer's direction, 6T 
is the spatial part of the Lagrangian temperature perturbation at 
the stellar surface and S ,, is the visible part of the stellar surface. 
The mode amplitude is normalized by 6To the root mean square 
of the perturbation over the whole stellar surface 



' g6THe,(p)dS] 



1/2 



(38) 



and the projected visible surface is normalized by nR^, the vis- 
ible disk surface for a star seen pole-on. With these normal- 
izations the disk-averaging factor of an hypothetical mode uni- 
formly distributed on the surface and seen pole-on is unity. The 
method of the calculation is detailed in Appen dix D and cor- 
rects t he calculation described in Appendix C of lLignieres et al.l 
(l2Q06h which actually only provides an approximate value of 
D(i). 

Figure |9] shows the spectrum of axisymmetric modes whose 
disk-averaging factor exceeds 2.5 percent. It appears that the 
disk-averaging effect does not allow to discard as many modes 
as for spherical stars. Indeed, in a given frequency interval and 
for the same visibility threshold, we find that the number of visi- 
ble modes is more than three time higher in the Q = 0.59Q.K star 
than in a spherical star Among the four classes of modes, the 
2-period island modes and the chaotic modes have similar visi- 
bilities and are significantly more visible than the 6-period island 
modes and whispering gallery modes. In Fig. |9] a few 6-period 
island modes are visible while no whispering gallery modes ex- 
ceed the chosen threshold. 

The relatively high visibility of the chaotic modes with re- 
spect to the 2-period island modes was not expected as the typi- 
cal horizontal wavelength of the chaotic modes is generally sig- 
nificantly shorter than the one of the 2-period island modes (see 
Fig. |6]l. We think that this is due to the iiTegular nature of the 
node pattern of the chaotic modes which makes the cancella- 
tion effect less effective than for regularly spaced nodes (like 
the whispering gallery modes). A practical consequence of this 
property is that, at such a rotation rate, methods to disentangle 
the regular spectrum from the iiTegular one should be developed. 
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Fig. 9. (Color online) Frequency spectra of axisymmetric modes 
where the amplitude is given by the normalized disk-averaging 
factor D{i) for a star seen pole-on / - and equator-on / = 
7t/2. Only frequencies such that D(i) > 2.5% are displayed 
and antisymmetric modes fully cancel out equator-on. The 2- 
period island modes (blue/dashed lines) and the chaotic modes 
(red/continuous lines) are the most visible while only a few 6- 
period island modes (magenta/dotted lines) and no whispering 
gallery mode exceed the 2.5% level. 



5.3. Seismic constraints 

Constraints on the star can be obtained once the island and 
chaotic modes sub-spectra are separated. Through the quanti- 
zation formulas of the regular modes, the asymptotic analysis 
provides the relation between regular frequency spacings and 
the physical properties of the star For example, according to 
Eq. (I3TI 1. the seismic observable 6„ depends on the value of the 
sound velocity along the path of the 2-period periodic orbit. The 
quantity 6( depends in addition on the second order transverse 
derivative of the sound velocity along the same path and the ra- 
dius of curvature of the bounding surfaces. The 2-period periodic 
orbit remains along the polar axis up to relatively high rotation 
rates (see Fig.[3]at Q/Q.K = [0.15,0.32] ). This implies that all 
the modes trapped within the corresponding island chain probe 
the center of the star which is known to be crucial for stellar 
evolution theory. It would be worth investigating if high order 
and low degree modes of rapidly rotating solar-type pulsators 
are in this case. Other informations on the star can be deduced 
from the numbers of island modes or chaotic modes. Indeed the 
number of mode in each class is related to the volume of the 
corresponding phase space regions (see subsection 4.7) which in 
turns depends on the stellar model. 

Further constraints on the star are expected from the identifi- 
cation of the chaotic modes. Contrary to the regular modes built 
on invariant torus, the modes built on chaotic region are not lo- 
calized in phase space and are expected be ergodic within their 
region of propagation. This property turns out to be of particu- 
lar interest for asteroseismology. The chaotic modes of the main 
chaotic region are indeed distributed all over the position space 
and do not avoid the stellar core like all the non-radial modes in 
non-rotating stars. Thus, in rapidly rotating stars, the chaotic p- 
modes have the potential to probe the physics of the core. While 
the sensitivity of the chaotic modes to this physics needs to be 
tested, quantum chaos studies indicate that the spatial distribu- 
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tion of chaotic modes is highly sensitive to ch anges in the models 
dSchack & Caveslll993HBenenti et al.Ll2002h . 

6. Conclusion 

We have constructed the ray dynamics in uniformly rotating 
polytropic stars and have presented the tools and concepts that 
enable to interpret it in terms of modes properties. Accordingly, 
the acoustic frequency spectrum of rapidly rotating polytropic 
stars is a superposition of frequency subsets associated with dy- 
namically independent phase space regions. The spectra associ- 
ated with near-integrable regions are regular while those associ- 
ated with chaotic regions are irregular but with specific statistical 
properties. Besides this global qualitative understanding of the 
spectrum organization, the ray dynamics also provides quanti- 
tative information. The EBK quantization of the near-integrable 
regions enables to explicitly construct the modes and the spec- 
trum from the ray dynamics. For the chaotic modes, the analysis 
of subsection 4.6 shows that there exists a parameter-free model 
of their frequency statistics. Moreover, we have seen in subsec- 
tion 4.7 that we can estimate from the ray dynamics the number 
of modes in each frequency subset. These results have been con- 
fronted with the properties of acoustic modes computed in the 
frequency interval [9wi, llcoi] showing that the present asymp- 
totic analysis does provide a quite accurate qualitative and quan- 
titative understanding of the actual spectrum in this frequency 
range. 

The analysis of section 5 argues for the importance of this 
asymptotic analysis for the mode identification and for the as- 
teroseismology of rapidly rotating stars. Indeed, the asymptotic 
results and the estimation of the mode visibilities tells us that 
the separation of the frequencies between chaotic and regular 
modes is a necessary prerequisite in order to analyze the spec- 
trum. When this is done, the observed regular spacings like 6„ 
and 6c can be related to the internal property of the star thanks 
to the asymptotic analysis. 

Further work on this theory could help the analysis of the 
observed spectra. First, it is important to establish more precise 
formulas such as Eq. ( |3TI ) for the regular modes corresponding 
to the different stability islands. The structure of chaotic modes 
at low frequency should be studied in more details, in frequency 
ranges lower than the ones used in the present work. This will 
allow to test the asymptotic analysis in frequency ranges where 
it is not supposed to hold in principle, but which are part of 
the observed spectra. Encouragingly, the regularity of the 2- 
period island modes has been already d emonstrated in a rela- 
tively low radial o rder range 5 < n < 10 jLignieres et al.Ll2006t 
iReese et al.L |2008|) . and more generally, in quantum mechanics 
the semi-classical analysis has been found to apply in energy 
ranges much lower than expected. Such a study would also allow 
to probe if the presence of scars (see section 4.6), which should 
be seen only at low frequencies, can organize part of the chaotic 
sub-spectrum. These studies will in particular enable to produce 
theoretical synthetic spectra which embody all the semi-classical 
information and can be used to test methods of analysis before 
dealing with actual spectra. 

Outside the scope of the asymptotic analysis per se, the 
mode identification would greatly benefit of accurate vis- 
ibility computations, of modes excitation studies and ob- 
viously of more re alistic models of centrifugally distorted 
rapidly rotating stars (Roxburgh', 120061: iMacGregor et all 120071: 
[Espinosa Lara & Rieutord, 2007). 

In conclusion, we believe that the asymptotic analysis we ex- 
posed is a promising way to interpret the spectrum of rapidly ro- 



tating stars. We have demonstrated that it can describe numerical 
spectra, and think that with suitable refinements it should pro- 
vide an important tool for the analysis of observed spectra such 
as those obtained by the instruments COROT and KEPLER. 
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Appendix A: The WKB approximation of the stellar 
oscillation equations ||6), 0, (H) 

The equations are first written in a normal form, then the eikonal 
equation is obtained using the WKB approximation. 



where 
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f( l-/)g() 

fsl 



/= 1 



The a function is given by: 
B 
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The expression of coc can be simplified in the limit w » A^o to 
give: 
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For polytropes Pq = Kg^^^^^'', the quantities describing the equi- 
librium model can be expressed in terms of the enthalpy ho as 
follows: 



go 



l^ho 



(A. 10) 



A.1. Normal form 



Equation (|A.9l l can then be simplified to give Eq. ( fTSl i. We note 
that while the w » A^o is expected to be valid for acoustic modes 
in real stars, the fact that the Brunt-Vaisala frequency becomes 
infinite at the surface of a polytropic model implies that this ap- 
proximation is not valid close to the surface of such models. 



We first eliminate the perturbation velocity u from the equations 
©, (|7|i, (O governing the evolution of the perturbations. Using 
equations (Q, the time derivative of Eqs. (|6]l and ([8]l reads: 



dip + V ■ (pgo) = AP, 



ciid^, + N^)p = di,P + ^go ■ VP 



(A.1) 



(A.2) 



A.2. Eikonal equation 

We look for wave-like solutions (|9]l of the normal form of the 
perturbation equation ( fTOb under the assumption that 1 /A the ra- 
tio between the wavelength of the solution and the background 
typical lengthscale is very small. Accordingly, the solution is ex- 
panded as: 

(5 = A(cl)„ + l(Bi..) A^Ao + {Ai.. (A.ll) 



where A^o is the Brunt-Vaisala frequency defined as 

Seeking harmonic solutions in time, the variable are written 
F = F expiiojt). Then, using Eq. ( IA.2b . p is expressed in terms 
of P and replaced in Eq. (lA.lb to yield: 



^2(l + £^)-c?/V.(^) 

2 J, 



P + 

go ■ VP - 



■ V(^o ■ VP) - oj'cifAP = 

so 



(A.4) 



We then look for a function a such that the substitution P - 
a*!* eliminates the first derivative term. The result is given by 
Eq. ([Tol l where : 



and the eikonal equation corresponds to the leading order of the 
expanded solution. 

The highest 0{A^) terms verify: 



2 2 



W2 
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+ ^A^(V(Do)l = A'iV^oY 



(A.12) 



where the (VOo)j_ - VOq - (VOq ■ mo)«o, «o being the outward 
unit vector in the direction opposite to the effective gravity. The 
effective eikonal equation then depends on the ordering of co/ojo 
with respect to A. If cij/loo - 0(A), then the above equation 
simplifies to: 



2 2 



= (VfDo)" 



(A. 13) 



which corresponds to the dispersion relation of acoustic waves. 
The ol>c term is retained because its sharp increase near the sur- 
face provokes the back reflection of the acoustic wave. 
On the other hand, if (l)/loo - 0(1), then we obtain 
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which corresponds to gravity waves. Note that this relation has 
been obtained under the assumption that the Coriolis force is 
negligible. While justified for high frequency acoustic waves, 
this assumption is not necessarily true for gravity waves whose 
frequency is limited by No- 

Finally, the next order of the expansion dA.l 111 yields th e am- 
plitude Ao in terms of "to (see for example lLighthilll (Il978l) ). 

Appendix B: Properties of the PSS 

Two specific properties of the rp{9) - rs(ff) - d PSS are con- 
sidered below. First, we check in the non-rotating case that the 
distance d can be chosen in order that all relevant trajectories 
intersect the rp(0) - rsiff) - d curve. Second, we define a coor- 
dinate system of the PSS which ensures that any surface of the 
PSS is conserved by the dynamics. 

6. 1 . Choice of the distance to the steiiar surface 

In the non-rotating case, equation (ISTT i enables to characterize 
the trajectories that do not cross the PSS for a given value of 
d, the distance to the stellar surface. The radius of the inter- 
nal caustic of these trajectories r,- must be such that r, > r^. 
Using the definition of r, and assuming that to » cjcird, im- 
plies that L > rpjcsifp). According to the relation between L 
and ( the degree of the coiTesponding mode L - (t + 1/2)/(l) 
(see Eq. (1231 ) in subsection 4.1), we find that these trajectories 
are associated with high degree modes > 136) for the cho- 
sen value of d - 0.08/?f and for w = 8.4wi, where wi is the 
lowest acoustic mode frequency of the stellar model considered. 
These modes are thus iiTelevant for asteroseismology since their 
amplitude strongly cancels out when integrated over the visible 
disk. 

B.2. Area-preserving coordinates of the PSS 

The PSS being defined by rp{6) - rsiO) - d, we show here that 
0, the colatitude, and kg = kg/cu), kg the angular component of k 
in the natural basis associated with the coordinate system - 
rs{6) - r, 6, 0], are area-preserving coordinates of the PSS. 

First we show that, for a general coordinate system x,-, the 
spatial coordinates x,- and the covariant component fe, of the 
vector k in the natural basis are conjugate variables of the 
Hamiltonian H given by Eq. ( fTTj i. The natural basis associ- 
ated to a coordinate system x, is defined by (Ei, E2, Ej,) where 
Ej = dx/dx'. The contravariant component k' of the velocity 
k = X thus verify x, = k' (the notation x,- denotes a deriva- 
tive with respect to the time-like coordinate f). If L(x, x, t) is the 
Lagrangian of a system expressed in a coordinate system x,, it is 
well known that a Legendre transformation enables to construct 
a Hamiltonian H - J] x/dL/dxi - L where - dL/dx/ is conju- 
gate to X,. The Lagrangian L - ^ - W{x) being associated with 
the Hamiltonian H given by Eq. the momentum variable pi 
can be simply computed: 

dL Idk k ~ dk ~ ^ ~ , , 

Pi = — = — =k—^kEi^ki (B.l) 

^ dxi 2 dk' dk' 

thus showing that [x,, ki\ are conjugate variables of H. 

Moreover, for a given conjugate coordinate system [x,, p,], 
the coordinates [x2, xj,,p2, P3] of the PSS de fined by xi = const 
are known to be area-preserving (Ott Jl993h . Thus, in our case, 
[^,0,k^,kg] is a conjugate coordinate system for the reduced 



Hamiltonian Hr and the system [0,kg] is area-preserving for 
the PSS ^ = const. Note that the natural basis and its con- 
jugate reads E^; = -er,Eg = (drs/dO) Cr + rpCg and E^ = 
-e,- + [{drJd0)/rp]eg,E^ = (l/rp)eg in terms of the unit vec- 
tor associated with the spherical coordinates {fir, eg). Thus, with 
respect to the wavevector components in spherical coordinates 
Tiy''^ ^lig'', the kg component reads kg - (drs/d0) fc*'''' + rpli^'' . 

Appendix C: Calculation of phase space volumes 

Following the Monte-Carlo quadrature method jPress et all 
1992), Ns points are randomly chosen in a known volume Vs 
that includes the volume V to be computed. The proportion of 
points inside V approximates the ratio V/Vs, thus providing an 
estimated value of V. The standard deviation error yields an es- 
timate of the relative eiTor, V(V'm/V - ^)/Ns, showing that the 
sampling volume Vs has to be as close as possible to the volume 
V and that the number of sampling points must be large. 

In our case, the main practical issue is thus to determine 
if a given point in phase space is inside or outside the 4- 
dimensional volume to be computed. The two volumes that 
we have computed are specified by two limiting frequencies 
9.42(^1 < H\.{x,k) < 11.85(1)1 and for each value of by the 
3D volume inside a given 2D torus. The first torus considered 
separates the whispering gallery region from the chaotic region, 
its imprint on the r - rp PSS being shown on Fig. lC.ll The vol- 
ume inside this torus includes the large chaotic region as well as 
the island chains around the 2-period and 6-period orbits. The 
second volume coiTesponds to the 2-period island chain and is 
delimited by a torus also shown on Fig. lC.ll 

To determine whether a given point jtq, is inside or outside 
these volumes, one could construct the r - rp PSS associated to 
the value of the Hamiltonian H^ixo, ko), advance the dynamics 
from xo, ko until the trajectory cross the r - rp PSS and find out 
whether the crossing point is inside or outside the torus. Here, to 
simplify the procedure, we used the fact that the r - PSS plot- 
ted against the scaled wavenumber k appeared to be insensitive 
to values of the frequency w in the domain considered. We thus 
consider the point Xo,ko/Hi.(X(),k()) and determine its location 
in the scaled phase space computed for a given frequency. To 
control this supposedly small frequency efi^ect the computation 
has been performed for the two extreme frequencies co - 9.42<yi 
and w = 11.85<x)i. Moreover, instead of advancing the dynam- 
ics up to the r - rp PSS, we construct a local PSS (either from 
a ( - const surface or a = const surface) to compare the 
location of the jco, ^0 point with the local imprint of the delimit- 
ing torus. In practice, the imprint of the delimiting torus is not a 
continuous curve as the torus is actually obtained from a space 
filling trajectory on the torus. We therefore follow such a trajec- 
tory over a sufficiently large number of time step to increase the 
number of point of the torus imprint on the different PSS. This 
procedure has been tested using trajectories which are known to 
be either inside or outside the torus (like for example trajecto- 
ries on nested tori inside the 2-period island chain). The number 
of points wrongly located by this procedure can be made small 
enough to have a negligible effect when compared to the statisti- 
cal error on the volume determination. Furthermore, the integra- 
tion domain has been divided in three sub-domains following the 
pseudo-radial direction ^. This enables to limit the ratio V/Vs as 
the sampling volumes can be more easily reduced in each sub- 
domain. 

Accordingly the number of modes within the 2-period island 
chain is 33 ± 1 if we use the bounding torus computed for oj = 
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Fig. C.l. (Color online) Intersection of two trajectories with the 
PSS at = 0.59f2/f . The crosses (blue) correspond to a trajec- 
tory on a torus containing most of the 2-period island chain. The 
diamonds (green) correspond to a trajectory on a torus bounding 
the central chaotic region. 



9.42(x)i and 34+1 for w - 1 l.SSwi. The effect of changing the 
frequency is small in this case and the number of modes can be 
estimated to 34 + 2. Likewise, the number of modes outside the 
whispering gallery region is 263 + 1 for uj - 9A2(±>\ bounding 
torus and 276 + 2 for (jj = 11.85^1. The frequency effect is larger 
but still reasonable for the present purpose. We took the mean 
value of 270 modes with an error of +8 mode. 



Appendix D: Calculation of the disk-integration 
factor 

According to the definition of the disk-integration factor, 
Eq. dJTI l. we are led to calculate integrals of the following form: 



/I. 
/I 



F{e, 0)dS ■ Ci 



G(0, 0, i)d^d^ 



(D.l) 
(D.2) 



where yu - cos and F{9, (f>) - W{0)e'""^ is the surface dis- 
tribution of the eigenfunction. The spherical coordinate system 
[r, 0, (p] is such that the polar axis is the rotation axis, and the 
meridional plane (p - Q \s chosen parallel to Ci. The condition 
dS • Ci = on the stellar surface defines a curve which separates 
the visible part of the surface S ,, from the invisible one (the sur- 
face is supposed to be convex). We call this curve the visibility 
curve. 

We used two methods to compute the integral /. The first 
one is approximate as it assumes that the visibility curve is con- 
tained in a plane but it is easy to compute accurately. The second 
method does not make this assumption but requires more com- 
puting time to complete accurate calculations. 

The vector dS at the star's surface reads: 



dS - sin Oks 



drs_ 

de 



2\ 



1/2 



e^d6d(f) 



(D.3) 



where denotes the unit vector perpendicular to the surface and 
rs{9) is the stellar surface. Thus, the function to integrate can be 
written as: 



where 



A ^ rs— (r,sm6)W(0) 
do 

B ^ r,^ if, cos 6) W{9) 
dO 



D.1. First method 



(D.5) 
(D.6) 



The colatitude verifying dS ■ Ci = for (p - Q is denoted 0l(O- 
As the inclination angle / can be restricted to [0, ;7r/2], we have 
7t/2 < Oiii) < n and the angle a defined as aii) - 6i{i) - 7r/2 
verifies < a{i) < n/l. We assume that the visibility curve is the 
intersection of the stellar surface with the plane sin ax -H cos az = 
0, where the vector Cq, of Cartesian coordinates (sin a, 0, cos a) 
is normal to this plane. 

Then, the integral is most simply calculated in the coordinate 
system in which the polar axis is aligned with the direction of 
the vector Ca- This coordinate system results from a rotation of 
angle a around the y axis of the original coordinate system, the 
new angular variables being denoted 0' and 0'. 

To express the integrand in these coordinates, we use the for- 
mula relating the spherical harmonics in both systems: 



+e 



(D.7) 



where d^^^,{a) do not generally have a simple form (Edmonds 
1960). Then, using the spherical harmonic expansion of G, we 
obtain: 



f=0 m=-C 



(D.8) 



(D.9) 



f=0 m=-C m'=-e 



Then, integrating over the longitude 0', from to 2ji, the terms 
involving Y'"'{9', 4>') vanish if m' + 0. It follows that 



e=o m=-f 

where we used the following relations. 



2(+l ' 



diid4> - du'difi' where ji' - cos 9' 
and defined Jc as. 



Je 



4jt 
2€+\ 



Y^(M')dM' 



if ( is even and ( i^O, 
if ^ = 0, 



(-l)^ 24"(m) if^isoddand^Ti 1, 
i if^=l. 



(D.IO) 

(D.ll) 
(D.12) 

(D.13) 
(D.14) 



To determine the coefficients G^„(i), we use the expression of 
G derived from Eq. ( ID.4l i: 



F(0,0)dS-ei = A(0) sine cos /e""*-B(0) sine sin/ cos (/.e'""^(D.4) G = A(9) cos ie'""^ - B(6») sin / cos 0e 



im(j) 



im(p 



(D.15) 
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It follows that 



D.3.1. Ellipsoid 



- \fki^m-\, m, m + 1 
so that the integral now reads: 



I/2n = I„-i + /,„ + /„,+! where 
= cos /A„,(q;) 
sin / ^ 

Im-\ —B,n-\{o:) 

sin / ^ 

hn+\ - — -;^B„,+\{a) 
where A,„ denotes: 



t=\m\ 



Al = In A(0)Y^(ff)smed0 



(D.16) 

(D.17) 
(D.18) 

(D.19) 
(D.20) 

(D.21) 
(D.22) 



the B„, terms being defined accordingly. 

Note that for modes which are equatorially anti-symmetric 
and axisymmetric (m - 0), Aoia) - JoA'^^Y^^ia) and Bi(a) - 

B-i{a) = 0, thus the integral / reduces to: 



I - An V^tAq cos(0. 



D.2. Second method 



(D.23) 



The visibility curve is no longer assumed to be planar. The in- 
tegration over the visible surface is first performed in the az- 
imuthal direction and then in latitude. If < < n/l - a , 
the integration is made between and 2n, while in the interval 
7r/2 - a < < 7t/2 + a one has to integrate between the two 
limiting azimuths -<Pl{0, and /) verifying dS ■ Cj = 0. The 
integration domain is thus divided in two sub-domains, such that 

5v. = [0, ^ - a] X [0,27r] u - a, ^ + a] X (D.24) 



According to Eq. (ID. 41 ). the integration over (p can be made 
analytically as it involves quadratures of e'""* and cos ^e'™** over 
[0, 27r] and {-4>l,iPl\- Then, depending on the value of m, the 
integral / reads: 



2n Jy" ad0 + 2 ^£_^ a(pi - b sm(cf)i)d0 if m = 

-n JJ'" bd0 + 2a sin(0z.) -b{(pL + ^ sin(20z.)) d0 
ifm-+l 

if TO ^ 0,±1. 



The surface being a quadric, the visibility curve is planar. 
Method 1 is therefore exact and should give the same result as 
method 2. In addition, the visible surface (obtained by taking 
F - I) can be obtained analytically. Indeed, the dimensionless 
equation of the ellipsoid is: 



Rp 

and in spherical coordinates: 
1 



1 + ^ C0S2 



(D.26) 



(D.27) 



where the distance have been normalized by the equatorial radius 
Re and Rp = Rp/Re- As dS is parallel to (2x, 2y, 2zIR],) and ei = 
(sin /, 0, cos /), surface points verifying dS ■ Cj = belong to the 
plane: 



. . zcos; 
xsini -\ — = 



By definition of the angle a, we have 



tan a - R^p tan / 



(D.28) 



(D.29) 



In addition, using Eq. ( ID. 281 ) and the relation between the 
Cartesian and spherical coordinates, the equation of the intersec- 
tion between the plane x sin/ +(z cos /)/^^ = and the ellispsoid 
is given by Eq. (ID. 27b and 



cot cot I 
coscpi = — -cot^cota 

R2 

p 

Thus, 

4>L{0,i) - arccos(-cot0cotci;) 



(D.30) 



(D.31) 



is defined if 7r/2 - a < < 7r/2 + a and verifies < (pi < n. 

The visible surface can be calculated analytically as the vis- 
ible curve is an ellipse. This can be seen using the Cartesian co- 
ordinates obtained by the rotation of angle a about the Oy axis: 



(D.32) 
(D.33) 
(D.34) 

The curve is then contained in the plane z' = and verifies: 



X — cos ax — sm az 

y = y 

z' — sin ax + cos az 



(cos a + 



(D.25) 



sin^ a 



-H = 1 



(D.35) 



^ sin[(/?i+l)0i] sm|(m-l 

m+l 



n[(m-l)fa] \^g 
m-1 / 

where a{0,i) - A{0) sin cos / and b{0, i) - B{0) sin sin ; . 



The surface of this ellipse can be calculated as well as its projec- 
tion in the direction Cj, denoted Sy : 



SiJR^ - 7TC0s(a - i) - 



1 + tan^ 



D.3. Tests 

The methods have been tested in the case of an uniformly dis- 
tributed function on the surface of an ellipsoid where they should 
both give the same result. Then, the error introduced by the ap- 
proximation of method 1 is estimated in the case of a Roche 
model surface. 



\ 1 + tan2 iR^ 
— 7TCOSI yjl + tan2 iRl 



(D.36) 
(D.37) 



Both method were successfully tested against this analytical 
expression. Method 1 is simpler as it does not require a numer- 
ical integration. It is also very accurate although, if one wants 
to reach machine precision, it is necessary to use the analytical 
value of 01. 
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D.3.2. Method 1 versus Method 2 

Method 1 is approximate as it assumes that the curve on the sur- 
face verifying dS ■ Ci = is planar and contained in the plane 
sin ax + cos az - 0. Nevertheless, the associated error is not ex- 
pected to be large as it concerns a small part of the whole visible 
surface. To test this, we simply compute the integral for F = I 
with both methods for a surface given by a Roche model of rotat- 
ing star. Method 2 computes the projected visible surface while 
the quantity calculated by method 1 is smaller as the integration 
includes a region where dS ■ Cj < and excludes a visible re- 
gion of equivalent surface which is symmetrical with respect to 
the star center The difference between the two quantities can be 
also calculated by directly computing the integral: 



a(f)L - b smi(pL)d0 (D.38) 



using either the correct value of (f)^ or the one corresponding to 
the assumption of method 1, that is 

cos (^'^ = -cot 6* cot a (D.39) 

In Fig. ID. II the relative error on the projected visible surface 
due to method 1 is plotted for Roche model surfaces of different 
flatness. 
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Fig. D.l. Relative error of the projected visible surface computed 
with method 1 for Roche models of different flatness : 0.1526 
(dot-dashed), 0.2594 (dotted) 0.2804 (dashed), 0.3092 (continu- 
ous line). 



It appears that, except for the near critical values of the flat- 
ness, the visible surface that is not considered by method 1 is a 
very small fraction of the total visible surface. Using method 1 is 
therefore a good approximation in these cases. For near critical 
flatness, the difference remains small although it can be useful 
to test results of method 1 with method 2. 



